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SELF-ADAPTIVE DIFFERENCE METHOD FOR THE EFFECTIVE SOLUTION OF 
COMPUTATIONALLY COMPLEX PROBLEMS OF BOUNDARY LAYER THEORY 

W. Schonauer, H.-G. Daubler, G. Glotz and J. Griming 
Flow mechanics research division 
Karlsruhe University computing center 

1. BACKGROUND 

As a spacecraft re-enters the earth's atmosphere (figure 1), /l * 

the air in the stagnation point area heats up so intensely that 
molecular oxygen and nitrogen dissociate and nitric oxide 
emerges as a significant intermediate compound. Consequently, 
chemical reactions take place in the air. A five-component air 
pattern would, for example, include the components N 2 , 0 2 , 

NO, N and 0. The composition is indicated by mass ruptures w^ 

= P^/p for all components i = 1 to N. 



Figure 1. Re-entry. 


♦Numbers in the margin indicate pagination in the foreign text. 
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The balanced equations; for total mass p, impulse 
(described by the velocity components u and v) , energy 
(described by temperature T) , and individual masses p^ are 
given to describe the flow. These balanced equations are 
supplemented by structural equations; for an example of energy 
current density vector or chemical production densities, see [1] . 

The resulting equational system is elliptical, i.e., a 
boundary value problem exists. For technically interesting 
cases, the calculation expenses for the numerical solution are 
insurmountable, even with the best computing machines. 

Therefore, an approximation theory for large Reynolds numbers is 
used instead, by which the flow is next calculated frictionless, 
and then transport operations are considered only in the 
boundary layer. In this way, heat transfer can be calculated on 
the contour. 


2. HYPERSONIC BOUNDARY LAYER EQUATIONS 

The resulting boundary layer equations and boundary 
conditions are given in [1], sections 7 and 8. Figure 2 shows 
the applicable coordinate system. As a result of the high 
flight altitude, the Reynolds number lies in the area at which, 
despite the high velocity, laminar flow exists. 



Figure 2. Boundary layer coordinates. 
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The non-linear differential equations do not have to be 
given here. The complete system is shortened to operator 
notation as follows: 

Equation of continuity % 

Impulse equation 
Energy equation 
Component equations of 
continuity for i=l to N > 

P x (operator P applied to x) indicates the entire equational /3 
system which is ordered so that the right side immediately 
becomes zero. The first component from P x indicates the 
equation of continuity, and so forth. One of the component 
equations of continuity is replaced by the relation 

N 

l w. - 1 = 0 (2) 

i=l 

since the system is otherwise linearly sloped. The system 
proves itself to be extremely "rigid" as a result of its 
inherent chemical reaction formula (small changes in the 
variables change the coefficients significantly) . The solution, 
with help from conventional differential procedures, is very 
"costly," because the necessary calculation of thermochemical 
functions for every grid point is extremely expensive to .compute. 


P X = 0 ; 


x = 


An 

v 

T 

w 


( 1 ) 


3. A FAMILY OF DIFFERENCE FORMULAS 

To solve the differential equations, we switch from the 
continuum to a difference grid (figure 3) ; derivatives are 
replaced by difference formulas, i.e. piece by piece the 
solution is represented by polynomials. The first derivatives 
appearing in the x-direction are replaced by a reversed 
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difference quotient, for example 


1 



Figure 3. Difference grid. 


This formula corresponds to the symbolism in the first line of 
figure 4: points indicate applicable function values and the 

circle indicates the place at which the difference quotient is 
formed. The symbol 0(h) means that a discretionary error arises 
from the order h. By adding further points, we get formulas of 
a higher order, which form a family of difference formulas. 

Since lengths that are not equidistant are accepted, we get, in 
the general case, an error of type O(h^), in which h indicates 
an "average" length. 

9 

' u x : * ' — & O(h) 'i. 04t»eR 

. — -±~® 0(h) 2, ORDER 

— * — * 0(h P ) F' le °*0 £*> 

Figure 4. Symbolism for difference formulas in x-direction. 
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It is important that the discrete error can easily be 
estimated by the application of a family of difference formulas 
with help from terms of a higher order of the family. If u xD 
or u _ are difference formulas of a higher or lower order of 
a family and or D ^ are the accompanying discrete errors, 
then the following applies, for example, to the derivative u : 

X 


u = u _ + D., = u _ + d_, from which follows 

x xD x 1 xD t 2' 


U xDj ^xD, 


(+ D„) , 


i.e., the discrete error can be estimated by the difference 
of the difference formulas up to an error . Yet this 
procedure is only significant as long as j | < D^j. If 
d 2 | > | Dj j, then the order of the difference formula 
(polynomial approximation) is "overdrawn." This "practical" 
estimation is of decisive importance for the controls of the 
operation, since the "exact" mistakes of type O(h^) are only 
calculable if the exact solution is known. The coefficients of 
the difference formulas, like the error formulas derived from 
them, are calculated by a sub-program with the help of 
interpolation polonomials directly at the computing machine. 


In the y~direction, the derivative u 
a family of central, equidistant difference 


is approximated by 
f ormulas ; . see figure 



at') 


Figure 5. a) central formulas, b) non-central formulas in 
the margin area. 
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5a. In the border region, asymmetrical formulas have to be 
(see diagram 5b) until the central formulas are applicable, 
most cases, we have the formulas of order ft. Derivatives 
u are likewise approximated; discretionary errors are 
estimated analogous to equation (4) by the difference from 
formulas of higher and lower order. 


4. DIFFERENCE METHODS 

The principle explaining how a non-linear system of 
differential equations is numerically solved should be 
illustrated by a simple model equation: 

P u = u - a (u) u = 0 , a > o . 

a yy 

Next, with the help of the Newton Process, a linearized 
"correction differential equation" is derived in which u is 
replaced by u v+1 = u v + Au v+1 , with v as the iteration 
numerator, which is nevertheless left out in the following 
equations. From (5) we get 

(u+Au) - a [u+ u] (u+Au) = 0 . 

A yy 

Through development from a to Au and linearization of the 
equation in Au, we get 

" Au x + a u U yy Au + a Au yy = u x ' a u yy t +0 < A u 2 )]. 

In operational notation, (6a) is shortened to 
Q Au = Pu , 

with Q as a linear differential equation operator. 


used 

In 


(5) 


(6a) 


(6b) 


Pu is indicated as the defect, which for an approximate 


/£ 
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solution u (=u V ) is not equal to zero. Equation (6b) 
illustrates a quadratic convergent iteration rule, by which /)u 
is to iterate until it has become "small enough." 

Equation (6b) is solved with the difference process, 
through which the derivatives appearing in Q and P corresponding 
to figures 4 and 5 are replaced by difference formulas. Figure 
6 shows the difference point which thereby arises and which 



Figure 6. Difference point. 

characterizes a dual parametric family of difference methods 
with the orders p and q as parameters. This implicit method was 
chosen in order to be able to control the length h^ 
independent of a stability condition. After the solution for 
the grid line x^ is calculated, we switch to line Xj, + ^, and 
so forth. The difference equation for the corrections Au^ on 
a grid line reads 

Qp. Au_ = (Pu)^ + D + D , (7) 

D D D x y 

with Q d as the matrix, which arises through the discretion of 

the operator Q and the related border region with (Pu) Q as the 

discretionary defect vector. D and D are the 

x y 

discretionary errors in the x- and y-directions . Solving from 
(7) to Au D (readily accomplished with the help of the 
Gaus-algorithm) yields 
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( 8 ) 




( p u) D 


+ D + D . 
x y 


We show equation (8) as an error equation. Up to the broken 
line, it represents the rule for the calculation of corrections 
for the Newton-Iteration, from which the defect (Pu) D is 
iterated to "zero." The subsequent calculated discretionary 
errors or become small if h or l becomes small. The 
solution error Au^, determined by the y-discretion, can be 
estimated by 


Au 


D;y 



D 


y 


(9) 


It is clear that, for effective computation, the three /7 

error terminals in the square brackets in (8) must be reasonably 
balanced. This will next be illustrated in general form. 


5. ERROR PROGRESSION 



Figure 7. Error progression. Key? 1) Re-entry problem; 
2) Accuracy demand; 3) Numerical result; 4) "Equations"; 
5) Numerical solution. 


In order to progress from the position of a technical ■ 
problem, for example heat transmission of a re-entering craft, 
to numerical results, we must run through many error-laden 


10 




single steps. The elements of the error progression are: 

1. Illustration of a mathematical model (using Navier-Stokes 
equations) . 

2. Simplification of the model (using boundary layer equations) . 

3. Coefficients from thermodynamics (for diffusion, reaction 
kinetics, etc.) 

4. Choice of a numerical solution method. 

5. Choice of calculation parameters (lengths). 

6 Discontinuance criteria for the iteration. 

The accuracy demand "presses" upon the error progression. 

The greatest error, that is the weakest member of the 
progression, determines the accuracy. It is not important to 
calculate with more numerical accuracy than that of the 
equational model. The calculation costs are minimized if alX 
numerical errors are of equal size and correspond to an 
established accuracy. 


6. TRANSFORMATIONS 

In the physical coordinate y, boundary layer thickness and 
dependent variables change tremendously with running length /8 
x; see figure 8a. By applying a modified Falkner-Skan 
transformation 



x 


/ u a (Od£ 



g = vp a 6 


( 10 ) 


we go from y, u, v, to y, f, g. As an exception, we might get 
"similar solutions," which no longer depend on x. In the usual 
case, the transformation reduces the change in the dependent 
variables considerably; see figure 8b. 
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Figure 8. Boundary layer a) 
transformation. 



without, b) with, similarity 


For the distant method, we are to divide the y-coordinate 
equally. Figure 9a shows a boundary layer profile, and figure 
9b adds to this the solution error, determined through the /£ 




Figure 9. a) Boundary layer profile? b) Error course; 
c) Asymptotic method. 


y-discretion. In the wall region, an error peak appears. By a 
suitable transformation we can extend the coordinate in the 
outer region (implicit grid angle) and reach the flattest 
possible error course. Figure 9c shows the influence of 
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coordinate y , from which boundary values are determined, on 
the boundary layer profile. A value for y that is too small 

cl 

produces large errors; a value for y that is too large wastes 
computation time. Therefore, the value of is chosen so 
that tg oi corresponds somewhat to the established accuracy of 
the solution. An implicit grid angle and normalization of the 
outer coordinates at the value one can be obtained by transition 
to the numeric coordinate y with the following transformation: 

y = Y a [cy + (l-c)y n ] , (11) 

whereby c: and n are so determined that the maximum error is 
minimized. Figure 10 shows transformation (11) and the grid 
angle in y direction which is then reached. 



Figure 10. Transformation equation (11). 


7. CONTROL IN y-DIRECTION 


The hypersonic boundary layer equations (1) now read in 
operator notation: 


P x = 0 




(u) 

(v) 

P i /P 


( 12 ) 
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As in equation (6) , we also give a Newton-iteration method 
(in this case at very great expense!) 

Q x = P x , (13) 

with the linear differential equation operator Q for the 
correction function Ax. Now we use the difference method /10 
with the difference point as in figure 6, and we get, as in (7) 
and (8) , the error equation 

lJ D - ^D 1 I(5?) D + *x + V ' (14) 

which is used to control the calculation. 

Figure 11 schematically shows the error determined by the 
y-discretion, for example temperature T for different orders q 
depending on length i. The dotted line shows the suitable 
combinations of X and q in each case. In order to find this 



Figure 11: Y-discret ionary error method for solution and 

determination of X and q. 
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combination for a given error F , the value of/, beginning 

9 

with q = 2, is determined, and this yields (l > L ■ ) as error 

F . Then we go to an order increased by about two (see arrows 

in figure 11) and so forth until the error grows as the order 

increases. Then the order would be "overdrawn" so that the last 

order is best. Figure 11 clearly shows that, for small errors 

F , a high order is best, but for large errors, a lower order 
9 

is better. This behavior very graphically answers the question 
"Is a high order better?" 


For boundary layer calculation, i-,q and y are determined 
by iteration of the beginning profile at the stagnation point. 
For the downstream calculation, JL and q remain fixed. This /II 
gives the error in (14) which is used to control the other 
errors . 


8. ITERATION DISCONTINUANCE AND CONTROL IN x-DIRECTION 


The defect (Px) D in (14) becomes quadratically smaller by 
successive iterations. An iteration step indicates calculation 
of a boundary layer profile. The iteration is then broken off 
as soon as 


IICPxU * 0.1 II D h i 

u Ko y Ko 


(15) 


exists, whereby Jl It indicates the maximum amount on a line x^ 
and Ko indicates a component-based comparison (per equation) . 
This makes the iteration error conform to D^,. 


In x-direction, length h. and order p have to become 
definite; see figure 6. Length h^ becomes so controlled on 
every point of x^ that 



(16) 
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takes effect. Thereby, D conforms to D . Order p is 

x y 

optimized by comparison with neighboring orders p+1: order p is 

optimal if the errors of the neighboring orders are larger 
(figure 12a) ; it can be increased if, for all components, the 
growing order error decreases (figure 12b) ; it must be lowered 
if, for one component, the error of order p-1 is smaller 
(polynomial "overdrawn"). 



Figure 12. Control of order p. 


In order to get good output profiles for iteration, the 
profiles for position are extrapolated with order p from 
the preceding profiles. The extrapolation represents /12 

a simple predictor. The complete solution method and its 
application to the hypersonic boundary layer equations are 
described in [2] - [5]. 



Figure 13. Extrapolation. 



9 . RESULTS 


Results are given for the flow on the so-called AGARD- 
hyperboloid, figure 14, namely for a typical case of unequal 
weight (A4) and a case of approximately equal weight (A) ; see 
[5]. The profiles for u/u^ as well as T/T a , in each case at 
the stagnation point, are given in figure 15. The profiles were 
calculated with 0.1% accuracy points for temperature T. The 
temperature profile in case A of approximately equal weight 



Figure 14. AGARD-Hyperboloid . 



Figure 15: Velocity and temperature profiles for AGARD- 

cases A and A4 at the stagnation point. 
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shows a "dent" caused by the chemical reactions. In figure .16, 
the profiles of mass rupture w^ are given at the stagnation 
point for the applied 5-component air pattern. 



Figure 16. Profile of mass rupture w^ for AGARD-case A4 
at the stagnation point. 

Figure 17 shows the progress of the Stanton number (heat 
transmission) and of the friction coefficient along the contour 
of the AGARD-hyperboloid for cases A and A4 (at 1% accuracy for 
temperature T) . 
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Figure 17a. Progress of Stanton number St along the 
contour coordinates x for cases A and A4 . 



Figure 17b. Progress of friction coefficient c^ along 
the contour coordinate x for cases A and A4. 


The computer program, produced with the support of the 
Volkswagen Foundation and the Society for Space Research, is 
available to those interested. It was modularly constructed so 
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that it can be easily altered for other flow problems or 
thermodynamic models. Detailed documentation is given in [6]. 

The calculations were done on the UNIVAC 1108 computer at the 
computing center of Karlsruhe University. 

10. NON- STATIONARY, 3-D, OR TURBULENT BOUNDARY LAYER 

The differential equations of the non-stationary hypersonic 
boundary layer read in operator notation as 

M "x , + Px = 0 (16) 

t 

with P and x from (12) and a diagonal matrix M. Inclusion of 
the additional time coordinate t means that, at the discretion 
in each moment tj , a complete two-dimensional boundary layer 
is to be evaluated; figure 18. 

/15 


Figure 18. Non-stationary boundary layer. 



In .the x- and t-directions there is a problem of initial value, 
and in the y-direction there is a boundary value problem. 
Accordingly, a completely implicit difference point is used 
(figure 19) with error orders p, q and r in the x-, y-, and t- 
directions. 
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Figure 19. 
layer . 


Difference point for non-stationary boundary 


The problem of generating an optimal difference grid is 
solved as follows: 

1. For the stagnation point x = 0, time t , from which the 
largest y-discretion error emerges, is determined by given 
values of JL and q. 

2. For t = t cr i t , stationary -£ opt / %pt determined at 
the stagnation point; then the stationary boundary layer 
along length x is solved, and from that lengths h^ are 
determined and accumulated. Then the x,y-grid is 
established. 

3. The non-stationary equations are solved, and from them time- 
length Tj and orders p and r are determined for each 
instant t^ analogous to section 8. Here the discretion 
errors are balanced so that 

M “l D xjM D y | < 17 > 

applies. 

Results from a multitude of examples show that, for /16 

chemically reacting boundary layers, the non-stationary terms 
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play a role only at very large time-gradients (as for the shock 
tubing method); see (7]. 

A self-adapting difference method for three-dimensional, 
incompressible boundary layer equations is being developed. 

Flow lines (x) and potential lines (y) are chosen as coordinates 
on the contour, and a similarity transformation is used. Error 
calculation shows that, by these means and with a relatively 
rough grid, good accuracy can already be achieved. The 
difference points are shown in figure 20 with optional error 
order in all three coordinates. 

Until now, these solution methods have only been used for the 
laminar case. Since the methods' are self-adaptive, they adapt 
easily to turbulent flow; in particular, by error calculation, 
they aid in finding grid angles which best capture the viscous 
substratum. But the problem of the turbulent boundary layer so 
far lies closer to element number 1 of the error progression 
from figure 7. As long as no reliable field equations for the 
turbulent flow are given, we cannot expect any exact results 
completely independent of applicable numbers. 




Figure 20. a) Half-implicit b) Fully-implicit , difference 
point for 3-D. 
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11. CLOSING COMMENTS 

Let's turn back one more time to the two-dimensional 
hypersonic-boundary layer and look closely at element number 3 
of the error progression in figure 7, which should contain the 
thermodynamic error. Calculation of the thermodynamic 
coefficients calls for a considerable portion of computing 
time. On one hand, the numbers do not have to be more exact 
than the accuracy made possible by thermodynamics. On the /17 
other hand, for some small given accuracy we could use a 
thermodynamic model that is simpler, less exact, and above all 
less expensive for making calculations. This problem gave rise 
to the following questions, which are currently being dealt with. 

1. How do the errors in the thermodynamic coefficients 
influence boundary layer results (for example, boundary 
layer profiles, heat transmission)? 

2. Which errors are produced by simple thermodynamic models, 
for example five reactions instead of eighteen reactions; 
air as a binary gas instead of a 5-component gas, etc.? 

Let us, in closing, still ask the question: Are such 

expensive procedures, as they are described here, really 
worthwhile? Next, as to computing time: the additional 

expenditure for error calculation and self-adaptation bring 
about a large-scale time savings. The expenditure pays off 
then, where very calculation- intensive problems are presented. 
Furthermore, because of error monitoring, we immediately get 
numerically reliable results without having to carry out 
expensive control calculations. Without self-adaptive methods, 
it would be nearly impossible to attain, for example, the exact 
calculation of non-stationary hypersonic boundary layers. Then, 
as for the programming expense: it is high, but here it is a 

matter of importance to program modularly. Many of our program 
building blocks (modules) of incompressible developments were 
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used for the chemically reacting boundary layer, the most 
modules of the stationary boundary layer being for the 
non-stationary , and so forth. If a versatile application is 
possible by skillful programming of flexible building blocks, 
i.e., a large part of the programming work is usable for other 
problems, then the high cost here is worthwhile. 
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